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A finite  difference  Lagranglan  code  calculation  of  the  firing  of  an 
earth  penetratcr  into  a soft  soil  site  has  been  carried  out.  The  penetrator 
was  a 400- lb,  6.9-ln.  diameter  steel  projectile  with  an  Initial  velocity  of 
900  ft/sec.  The  field  tests  were  conducted  In  July  1974  at  the  DRES  Site  in 
Canada.  The  penetrator  was  treated  as  a rigid  body  In  the  code  solution.  A 
four-layer  idealization  of  the  target  site  was  used.  The  calculated  decel- 
erations were  approximately  110  g's,  I29  g's,  and  49  - 90  g's  in  the  first 
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20.  ABSTRACT  (Continued) 


threu  luyors.  Frictional  forces  at  the  peneti'ator/soil  interface  were 
responsible  for  a major  part  of  these  decelerations.  (The  calculation 
was  not  carried  into  the  fourth  layer.)  Eiqjerimentally-observod  decel- 
erations were  somewhat  smaller  than  the  calculated  values,  but  the  general 
agreement  between  calculation  and  experiment  was  good.  It  is  concluded 
that  finite  difference  code  solutions  can  provide  useful  information  about 
earth  penetrators,  particularly  for  those  portions  of  penetration  events 
where  loads  and  decelerations  are  changing  rapidly. 


PREFACE 


This  report  describes  a finite  difference  code  calcu- 
lation of  an  earth  ponclrator  experiment  conducted  in  July 
1974  at  the  DRES  site  in  Canada.  The  work  was  carried  out 
within  a short  time  schedule  in  an  effort  to  provide  a pre- 
shot calculation.  The  principal  investigators  at  : 
California  Research  9 Technology  were  M.  H.  Wagner  and 

K.  N.  Kreyonhagen.  W.  S.  Goorke  performed  the  computer 
programming  and  code  development  work.  C.  C.  Pulton  and 

L.  A.  DeAngelo  provided  valuable  assistance  in  formulation 
of  the  material  models  and  data  reduction,  and  in  the  com- 
puter production  runs,  respectively. 

The  site  profile  and  material  properties  were  provided 
by  J,  Cl.  Jackson,  Jr.,  1).  S.  Army  Engineer,  Waterway. s 
Experiment  Station. 

The  Project  Officer  was  Major  T.  1).  Stong,  SPSS, 
Uol'en.sc  Nuclear  Agency. 
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SECTION  I 


INTRODUCTION  AND  SUMMARY 

1.1  BACKGROUND  AND  APPROACH 

For  some  applications  of  both  nuclear  and  high  explo- 
sive systems,  it  i.s  desirable  to  have  the  weapon  penetrate 
into  an  earth  target  before  exploding.  Ample  demonstra- 
tions have  been  made  which  show  that  substantial  penetra- 
tions are  achievable,  at  least  in  soft  media,  using  high 
velocity  penetrators  (o.g.,  Ref.  1). 

The  objective  of  the  present  study  was  to  examine  the 
current  capabilities  of  finite  difference  codes  to  predict 
and  analyze  earth  penetrator  events.  This  was  done  by  cal- 
culating a specific  penetration  event,  which  will  be  re- 
ferred to  herein  as  the  DNA  earth  penetrator  experiment. 

The  calculation  was  performed  during  June  and  July  1974. 
Specific  information  sought  from  the  calculation  was  the 
deceleration,  velocity,  and  displacement  histories  for  the 
penetrator,  the  surface  loading  history  on  the  penetrator, 
and  a description  of  the  earth  target  response. 

The  DNA  earth  penetrator  experiment  was  performed  in 
July  1974  by  Sandia  Laboratories,  Albuquerque.*  This  ex- 
periment consisted  of  a series  of  field  te.sts  conducted  at 
the  Watching  Hill  Blast  Range  of  Defense  Research  Estab- 
lishment Suffield  (DRES),  Alberta,  Canada.  The  tests  in- 
volved vertical  firings  of  instrumented  penetrators  at  500 
ft/sec  into  the  ground.  The  penetrators  were  steel  cylin- 
ders with  ogival  noses,  weighing  about  400  lbs,  with  a 
total  length  of  about  five  feet  and  a diameter  of  6.5 
inches, 


The  curth  target  metUa  at  the  test  site  is  compara- 
tively uniform  to  the  depths  of  interest,  consisting  of 
layered  glacial  lake  deposits.  Using  extensive  soil  test 
data  previously  gathered  in  the  general  area,’  recommended 
properties  for  a four-layer  idealized  site  profile  were 
prepared  by  the  Waterways  lixperimont  Station  (WliS)  . This 
information  was  used  to  develop  tho  material  models  for 
the  code  calculation.  Tho  sketch  below  indicates  the 
general  conditions  for  modeling  the  DNA  earth  penetrator 
e)^periment. 
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Figure  1.  Calculational  Model  of  DNA 
Earth  Penetrator  Experiment 
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The  calculation  described  in  this  report  was  made 
with  the  WAVE-L  code,  a two-dimensional  Lagrangian  finite 
difference  technique  based  on  the  HEMP  method,''  The  pene- 
trator  was  represented  as  a rigid  body,  which  is  a realis- 
tic approximation  when  the  penetrator  deformations  are  en- 
tirely elastic  . 

In  addition  to  the  code  analysis  performed  by  CRT 
which  is  described  herein,  separate  calculations  of  the 
same  event  were  performed  by  Sandia  Laboratories, 
Albuquerque*  and  by  Physics  International.* 

1.2  SUMMARY  OF  RESULTS 

Figure  2 shows  the  calculated  axial  force  and  deceler- 
ation of  the  penetrator  vs  depth.  The  penetrator  problem 
was  integrated  completely  through  Layer  1.  After  an  ini- 
tial buildup  (as  the  nose  was  entering  the  soil),  quasi- 
steady state  conditions  were  established  and  the  penetrator 
deceleration  remained  fairly  constant  at  about  110  g's. 

This  increased  to  about  125  g’s  after  the  nose  had  fully 
entered  Layer  2. 

Figure  3 shows  the  Lagrangian  grid  with  the  penetrator 
tip  at  200  cm  depth  (6.6  ft),  when  it  was  completely  in 
Layer  1.  Soil  separates  from  the  penetrator  about  15  inches 
back  from  the  nose  tip.  Largo  frictional  forces  along  the 
penetrator-soil  surface  cause  severe  distortion  of  the 
cells  adjacent  to  the  penetrator.  Periodic  local  rezonlng 
of  the  grid  was  required  in  this  region. 

To  avoid  the  time-consuming  and  largely  uninformative 
integration  through  Layer  2,  the  solution  was  stopped  after 
the  nose  had  fully  entered  Layer  2,  and  restarted  when  the 
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Figure  2.  Axial  Force  on  Penetrator  (showing  total  axial  force  and 

frictional  component)  and  Deceleration  of  Penetrator  vs  Depth 
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nose  reuchcJ  Luyer  3.  The  penotrator  velocity  was 
reduced  to  reflect  a constant  125-g  deceleration  through 
Layer  2. 

As  the  penetrator  entered  Layer  3,  the  deceleration 
dropped  sharply,  reaching  4S-50  g's  after  the  nose  was 
fully  embedded  in  Layer  3.  Assuming  the  deceleration 
would  remain  constant  at  47  g's  in  Layer  3,  it  was  esti- 
mated that  the  penetrator  would  reach  Layer  4 (at  24  ft 
depth)  at  59.5  msect  and  with  a velocity  of  330  ft/sec, 
or  101  m/sec.  The  solution  was  not  extended  into  Layer  4 
due  to  time  constraints.  . 

As  is  seen  in  Figure  2,  friction  was  responsible  for 
approximately  half  of  the  total  deceleration  in  the  calcula- 
tion in  all  three  layers. 

Figure  4 shows  a comparison  of  the  calculated  decel- 
eration history  with  preliminary  experimental  results. 

The  large  excursions  and  the  absence  of  obvious  demarca- 
tions between  layers  in  the  experimental  results  once 
again  confirm  the  difficulty  of  using  Idealized  soil  pro- 
files to  represent  actual  geologic  conditions.  The  some- 
what high  i.evel  of  deceleration  predicted  by  the  calcula- 
tion in  each  layer  suggests  that  excessive  friction  was 
specified  in  the  model.  Despite  these  problems,  the  agree- 
ment shown  in  Figure  4 is  regarded  as  very  good  for  an 
initial  demonstration  problem.. 

The  basic  zoning  for  the  calculations  used  4 cells 
per  penetrator  diameter.  The  tcTr.D..ter  time  (exclusive  of 
bad  runs,  check  problems,  and  ed‘.\.s}  was  about  12  hours  on 
the  Univac  1108. 
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1.3  CONCLUSIONS 


The  Lagrangian  finite  difference  calculation  of  the 
DNA  penetrator  problem  produced  generally  good  agreement 
with  experimental  acceleration  records.  This  indicates 
that  the  WAVE-L  code  can  provide  useful  information  re- 
garding important  aspects  of  rigid-body  penetration  into 
earth  targets.  Results  obtained  from  another  program’  in- 
dicate that  WAVE-L  can  also  be  used  to  determine  stresses 
and  accelerations  within  earth  penetrators.  Thus,  finite 
difference  codes  can  be  applied  to: 

o Studying  phenomenology  of  penetration  in  differ- 
ent earth  media  types 

o Examining  the  effects  of  variations  in  target  pro- 
perty, penotrator  design,  and  impact  parameters 
upon  the  penetration  processes 

0 Establishing  the  design  environment  (stresses 
and  accelerations)  within  an  earth  penetrator 
in  representative  types  of  targets  and  under 
different  impact  conditions. 

Problems  undeniably  exist  in  the  application  of  finite 
difference  codes  to  earth  penetrator  analyses.  However, 
with  reasonable  assurance,  some  of  these  can  be  overcome 
with  modest  additional  effort.  The  others,  we  argue,  can 
be  avoided  by  recognizing  the  objective  of  the  numerical 
solutions . 

In  the  former  category,  the  needs  for  realistic  char- 
acterization of  friction,  for  efficient  correction  of  grid 
distortion,  and  for  establishment  of  zoning  criteria  for 
reasonably  accurate  solutions  should  all  yield  to  current- 
ly planned  work. 
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Tho  more  difficult  problems  have  to  do  with  long 
computer  running  times  and  with  persistent  geologic 
uncertainties . 

Long  running  times  are  involved  in  attempting  to  solve 
complete  deep  penetration  problems,  and/or  in  predicting 
depths  o£  penetration  in  soft  targets.  In  such  problems, 
however,  major  portions  of  the  penetration  may  involve 
phe^nomena  which  change  only  slowly  with  time.  It  is  not 
necessary  or  fruitful  to  analyse  such  quasi-steady  state 
problems  with  a finite  difference  code.  Rather,  such  codes 
should  be  reserved  for  analyzing  those  portions  of  penetra- 
tion events  where  loading  on  the  penetrator  is  most  severe 
or  is  changing  rapidly  (e.g.,  while  the  penetrator  nose  is 
entering  the  surface,  or  whilu  it  is  passing  between  dif- 
ferent geologic  layers,  or  perhaps  just  before  it.  comes  to 
rest) . 

Some  geologic  uncertainties  will  always  exist,  even  when 
extensive  svibsurface  surveys  are  possible.  However,  such 
uncertainties  are  generally  troublesome  only  when  one  is 
attempting  to  validate  a predictive  method  by  comparing  it 
with  a test  at  a specific  site.  Exact  knowledge  of  real 
target  properties  is  unattainable,  nor  is  such  knowledge 
needed  for  the  primary  code  applications  listed  at  the  be- 
ginning of  this  section. 


SECTION  II 

DNA  PENBTRATOR  PROBLEM  DESCRIPTION 

The  overall  problem  conditions  for  the  numerical  solu- 
tion are  shown  in  Figure  S. 

2.1  PENBTRATOR 

The  penetrator  was  a steel  cylinder,  with  an  ogival 
nose  section  and  a small  45^  conical  tip,  as  Indicated  In 
the  insert  in  Figure  S.  For  this  code  solution,  the  pene* 
trator  was  treated  as  a rigid  body.  This  assumption  is 
reasonable  for  analyses  of  penetration  dynamics,  inasmuch  as 
the  penetrator  undergoes  only  small,  elastic  deformations 
which  do  not  significantly  interact  with  the  penetration 
processes  in  the  target.  Separate  analyses,  either  finite 
difference  or  finite  element,  can  be  performed  to  deter- 
mine the  stress  environment  within  the  projectile,  as  is 
illustrated  in  Section  VI  of  this  report. 

2.2  TARGET  MEDIA 

The  target  media  were  modeled  to  correspond  with  the 
Watching  Hill  site  at  DRES,  Alberta,  Canada.  The  soil  pro- 
file was  idealized  by  WES  into  four  distinct  layers.  Data 
from  soil  property  measurements  by  WES  were  used  as  the 
basis  for  constructing  the  models  for  each  layer. 

The  basic  properties  used  for  the  four  layers  are 
listed  in  Table  1.  Comparative  uniaxial  strain  load- 
unload  paths  for  each  layer  are  shown  in  Figure  6.  The 
failure  surfaces  for  each  layer  are  compared  in  Figure  7. 


Layar  1.  P.  ■ 1.4897  gm/cm  ■ 9 


Layer  2.  ■ 1.42S6  gm/cm"  ••  89  pcf 


Layer  3.  « 1.8381  gm/cm  11|6  pcf 


Layer  4.  p • 1.9702  gm/cm"’  • I2j3  pcf 


Figure  5.  Problem  Condltlona  for  CRT  Rigid  Body 

Solution  of  DMA  Earth  Penetrator  Experiment 
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Layer  4 is  considered  to  be  water  saturated,  with  a water 
content  of  '301,  accounting  for  its  low  shear  strength. 

For  all  layers,  a hystoretic,  elastic- ideally  plastic 
model  with  a non-associatod  flow  rule  (Prandtl -Reuss)  was  I 

formulated.  In  loading,  the  bulk  and  shear  moduli  are  | 

functions  of  volumetric  strain  (y  • p/Pq) • unloading,  5 

the  hulk  modulus  is  a function  of  maximum  volumetric  strain 
(Umgjj)  or  maximum  pressure  The  unloading  shear 

modulus  is  constant  in  Layers  1,  2,  and  3.  In  Layer  4,  the 
unloading  shear  modulus  Is  dependent  on  the  unloading  bulk 
modulus.  Variable  failure  surfaces,  where  the  yield  condi- 
tion is  a function  of  the  pressure,  were  employed. 

Detailed  descriptions  of  the  material  models  are  con- 
tained in  Appendix  A. 

2.3  FRICTION  ALONG  PBNETRATOR- TARGET  SOIL  INTERFACE 

I A general  friction  rule  was  formulated  and  programmed 

into  the  code  that  includes  frictional  stresses  In  the 
equations  of  motions  of  both  sides  of  the  sliding  Inter- 
. face.  The  code  can  readily  accept  any  friction  rule. 

For  this  problem,  a simple  friction  rule  was  used  to  com- 
. puto  the  applied  shear  stress  on  tho  soil-penotrator  inter- 

face; l.c. , 

T ■ . 

I where  Y Is  the  yield  surface  equation  of  the  medium  and 

\ Is  tho  stress  component  normal  to  tho  ponetrator  sur- 

I face.  The  yield  surface  equations  for  each  layer  are 

I given  in  Appendix  A.  Selection  of  the  .6  coefficient  was 
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based  partially  on  discussions  with  WES.  This  friction 
rule  produced  relatively  large  frictional  stresses,  which 
accounted  for  40*65%  of  the  total  decelerating  force  in 
the  first  three  layers. 

Time  constraints  did  not  permit  the  use  of  a more 
sophisticated  friction  model  for  this  calculation.  Addi* 
tional  work  is  underway  which  accounts  for  limitations  due 
to  the  internal  friction  in  the  soil  and  the  effects  of 
failure.  Guidance  from  experiments  measuring  frictional 
parameters  under  conditions  which  are  relevant  to  earth 
penetrators  is  also  needed. 
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SECTION  III 


COMPUTATIONAL  METHOD 

The  WAVE-L  code  was  employed  for  this  calculation. 
WAVE-L  is  a two-dimensional  code  which  solves  the  equa- 
tions of  motions  for  elastic-plastic  bodies  by  means  of  a 
finite  difference  Lagrangian  technique.  The  mathematical 
formulation  is  basically  the  same  as  that  described  by 
Wilkins. ^ Lagrangian  cells  are  normally  attached  at  the 
lattice  points.  In  order  to  allow  penetration  of  the  pro- 
jectile through  the  computational  grid)  a decoupling  cri- 
teria was  used  along  the  axis  of  symmetry  in  the  soil  med- 
ium. A sliding  interface,  in  conjunction  with  n friction 
model,  was  used  between  the  rigid  body  and  the  soil.  A 
lithosttttlc  field  in  the  soil  was  generated  by  proscribing 
initially  compressed  soil  states. 

3.1  BASIC  GRID  DESIGN 

The  rigid  contour  of  the  penetrator  periphery,  includ- 
ing the  conical  tip,  was  resolved  with  a total  of  44  lat- 
tice points.  The  spacing  in  the  vertical  direction  began 
with  .5  in.  at  the  tip,  gradually  Increasing  to  1.625  in. 
along  the  ogival  .section,  and  continuing  at  a uniform 
1.625  in.  to  the  rear  of  the  projectile. 

The  computational  grid  in  the  soil  was  divided  into 
two  principal  regions,  as  shown  in  Figure  8.  From  the 
axis  of  symmetry  out  to  19.5  in.  radius,  all  cells  were 
initially  uniform,  with  ■ 1.625  in.  (Thus  there 

are  Initially  4 cells  per  penetrator  diameter.)  Beyond 
19.5  in.  radius,  Ar^  increases  in  lOt  steps.  The  grid 
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extended  to  a radius  of  7.1  ft  (bo/ond  the  field  of  view 
in  the  figure).  In  this  outer  region,  also  increases 

(with  increasing  depth)  in  101  steps  until  i“  reaches  a 
maximum  value  of  . This  design  for  the  LUter  region 

conserves  computational  cells,  without  significantly  de- 
grading the  solution. 

A slide  lino  was  specified  that  runs  along  the  axis 
in  the  soil  and  along  the  penetrator  contour.  As  the  pene- 
trator  impinges  into  the  target,  this  dual  line  defines 
the  penetrator- soil  interface  and  allows  slippage  and/or 
separation  to  occur.  Another  slide  line  is  required  along 
the  boundary  between  the  two  grid  regions.  This  line  is 
needed  to  accommodate  the  differing  alignment  of  the  grids 
on  both  sides,  but  no  substantial  slippage  occurs. 


3.2  "ROLLING  GRID"' 


If  the  initial  computational  grid  is  required  to  cover 
the  entire  target  region  which  is  involved  in  the  complete 
penetration  event,  the  grid  will  involve  an  extraordinary 
numer  of  cells  in  order  to  achieve  reasonable  spatial  re- 
solution. The  alternative  a smaller  number  of  large 
cells  (if  the  coarse  resolution  is  acceptable).  In  either 
event,  the  efficiency  of  a solution  of  an  event  (such  as  a 
deep  penetration)  where  the  action  is  predominantly  in  one 
direction  can  be  substantially  Increased  by  deleting  cells 
after  the  "action"  has  passed,  and  adding  new  cells  in  the 
region  where  f»ctivity  has  not  yet  started.  This  was  done 
in  the  current  solution  by  deleting  rows  of  cells  in  the 
target  material  behind  the  penetrator  (where  the  material 


was  judged  to  have  little  or  no  effect  on  the  continuing  sol- 
ution) and  simultaneously  adding  rows  of  cells  in  the  as-yet 


inactive  region  ahead  of  the  penetrator.  This  technique 
conserve*:  computational  colls  and  thus  increases  the 
available  spatial  resolution  which  is  attainable  with  a 
given  number  of  cells  in  a problem.  Thirty-two  vertical 
grid  lines  and  a maximum  of  106  horizontal  grid  lines  were 
used.  The  number  of  computational  cells  in  the  problem  ;<t 
any  one  time  ranged  from  1300  to  1800. 

3.3  PRE-DEPORMED  GRID  AND  LOCAL  REZONING 

During  penetration,  the  cells  in  the  soil  near  the 
penetrator  become  severely  distorted.  This  distortion 
occurs  for  two  main  reasons:  (1)  cells  from  near  the  axis 

are  squeezed  as  they  are  pushed  outward  in  the  flow  around 
the  penetrator  nose,  and  (25  where  large  frictional 
stresses  are  present  along  the  sliding  interface,  cells  are 
distorted  by  the  shear.  While  such  distortion  appears  physi- 
cally cdnsi-steut  with  the  frictional  rule  ann  the  penetration 
dynamics,  it  causes  computat  tonal  inaccura  .ic.s  and  reduces 
the  integration  time  step. 

To  correct  extreme  distortion,  local  rezones  were 
periodically  performed  during  the  course  of  the  integra- 
tion. In  a local  rezone,  the  computational  grid  in  the 
distorted  region  is  repositioned  so  as  to  give  met  regular 
cell  shapes.  A comprehensive  rezone  processor  i.*.  used  for 
this  purpose  that  correctly  redistributes  the  cull  varia- 
bles among  the  new  cells. 

To  increase  the  time  interval  between  loca.’.  rezones, 
a special  grid  was  designed,  in  wit  .ch  cells  in  the  f irst 
two  columns  near  the  axis  are  canted  upward.  As  the  pone- 
trator  reaches  these  cells,  fr:.ctional  forces  deform  the 
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cells  downward,  The  pre-derormation  is  a geometric  con- 
venience only,  and  does  not  affect  the  stress  field  which 
is  generated  in  any  important  way.  Figures  9 and  10  show 
the  grid  before  and  after  a local  rezone.  Note  in  Figure 
9 that  the  pre-deformed  cells  ahead  of  the  penetrator  be- 
come "bent  down"  as  they  drag  along  the  penetrator  surface. 

3.4  GRID  DECOUPLING 

The  sliding  interface  formulation  in  WAVB-L  includes 
the  capability  for  dynamically  decoupling  or  unlocking  the 
grid  points  as  a prescribed  decoupling  or  failure  criter- 
ion i^  met.  Decoupled  grid  points  can  slide  along  the  in- 
terface and  also  separate  (forming  a void).  In  a penetra- 
tion problem,  decoupling  permits  the  cells  initially 
attached  to  the  axis  to  flow  around  the  advancing  nose 
shape.  In  this  calculation,  decoupling  occurred  when  the 
generalized  plastic  strain  reached  It.  Generalized  plastic 
strain  is  a measure  of  the  cumulative  plastic  deformation 
an  element  undergoes,  and  is  defined*  as 

e P ■ / . 

_p  «.p 

de  is  always  >0  , thus  e monotonically  increases  as 
any  type  of  plastic  strain  is  experienced. 

3.5  VOID  OPENING  AND  CLOSING 

During  high-velocity  penetration,  the  soil  flows 
around  the  nose  and  typically  separates  from  the  penetra- 
tor before  reaching  the  tangenjy  point.  This  separation 
and  the  resultant  void  opening  are  critical  factors  in  the 
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Figure  10.  Grid  Configuration  after  Local  Rezone 
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calculations,  since  they  determine  the  effective  contact 
area  over  which  frictional  forces  are  applied.  In  this 
calculation,  separation  of  a soil  point  could  occur  when 
the  stress  normal  to  the  penetrator  was  no  longer  compres- 
sive; i.e,,  ^ ^ ' When  this  criterion  is  met,  the 

point  is  moved  as  a free  surface  point.  (The  soil  point 
may  subsequently  collide  with  the  penetrator,  in  which 
case  it  again  becomes  attached  to  the  penetrator.  Closing 
of  the  void  presumably  occurs  during  the  terminal  stages 
of  penetration.) 

3.6  "LEAPPROOOlNa" 

During  the  penetration  of  Layer  1,  it  was  observed 
that  the  deceleration  remained  at  a relatively  constant 
level  following  the  burial  of  the  penetrator  nose.  This 
should  be  the  case,  so  long  as  the  penetrator  velocity  Is 
not  drastically  changed  and  if  the  contact  area  of  the 
nose  with  the  soil  remains  about  the  same.  Since  no  im- 
portant new  information  la  gained  by  integrating  out  the 
complete  penetration  of  a layer,  it  was  decided  to  skip 
ahead  after  the  penetrator  nose  had  become  fully  Immersed 
in  a layer.  To  accomplish  this,  the  solution  is  stopped 
after  the  deceleration  level  in  a layer  has  been  estab- 
lished. It  is  restarted  (with  an  appropriately  reduced 
velocity)  with  the  nose  a short  distance  from  the  next 
layer.  This  technique  is  illustrated  in  Figure  11. 
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Aicar  noit  tuiiy  DuriMi 
dtcaltratlon  In  Layer  1 
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'Phaaa  I tamlnatei  after 
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Figure  11.  Illustration  of  "Leapfrog"  Technique  for  Computing  Deep  Penetration 


SECTION  IV 


NUMERICAL  RESULTS  AND  EXPERIMENTAL  COMPARISONS 
4.1  PENETRATOR  DYNAMICS 

The  penetration  was  nuroerically  integrated  completely 
through  Layer  1,  and  into  Layer  2 to  a depth  of  300  cm,  or 
9.8  ft.  (This  la  identified  as  Phase  I in  Figure  11.) 
Figure  12  shows  the  computational  grid  at  the  end  of  Phase 
I.  (Badly  distorted  cells  have  been  deleted  from  the  in- 
side of  the  hole  behind  the  penetrator.  These  cells  no 
longer  affect  the  penetration,  and  their  removal  improves 
computational  efficiency.) 

Figure  13  shows  the  axial  force  on  the  penetrator  and 
the  penetrator  deceleration  vs  depth.  The  deceleration 
climbed  sharply  as  the  penetrator  nose  first  entered  the 
soil,  reaching  a level  of  about  110  g's  where  it  remained 
fairly  constant  during  the  remainder  of  Layer  1.  As  the 
penetrator  entered  Layer  2,  the  deceleration  Increased 
somewhat,  reaching  about  125  g's.  Based  on  the  experience 
in  Layer  1,  it  was  assumed  that  the  deceleration  would  con- 
tinue at  125  g's  in  Layer  2.  This  assumption  forms  the 
basis  for  the  "leapfrogging"  approach  described  in  Section 
3.S,  by  which  the  time-consuming  and  largely  uninteresting 
integration  of  the  quasi- steady  state  penetration  within  a 
geologic  layer  is  bypassed.  Phase  1 of  the  calculation  was 
thus  stopped  with  the  penetrator  at  300  cm  (9.8  ft)  depth. 
This  is  56  cm  into  Layer  2.  Its  velocity,  , at  this 
time  was  131  m/sec  (430  ft/sec). 
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"x"  In  call  means  material  is  currently  on  yield  aurface  with  P>0; 
means  material  falling  with  P<0. 


Figure  12. 


Orld  Configuration  at  End  of  Phase  I of 
the  Solution,  Penetrator  at  300  cm  Depth 


. .1  i..‘  l..±' 
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Flgura  13.  Axial  Forca  on  Panatrator  (ahowing  total  axial  force  and 

frictional  conponant)  and  Dacalaratlon  of  Panatrator  va  Depth 


Deceleration  of  Penetrator 


To  start  Phase  II  of  the  calculation,  the  penetrator 
depth  was  set  at  469  cm,  or  15.4  ft,  and  it  was  given  an 

initial  velocity,  , of  114  m/sec,  or  374  ft/sec. 

(This  was  determined  by  applying  12S-g  deceleration  to 

, the  final  velocity  at  the  end  of  Phase  I.)  The  ini- 
tial particle  velocities  in  the  soil  field  around  the  pene- 

trator  for  Phase  11  wore  set  by  reducing  the  field  veloci- 
ties at  the  end  of  Phase  1 by  the  ratio  Vjj  /Vj^  ■ .87  . 
This  gives  a reasonable  estimate  of  the  field  velocities, 
and  inaccuracies  are  quickly  corrected  as  the  integration 
continues. 

Using  these  starting  conditions,  the  Phase  II  inte- 
gration was  carried  out  as  the  penetrator  nose  moved  from 
Layer  2 into  Layer  3,  to  a depth  of  536  cm,  or  17.6  ft. 

As  is  seen  in  Figure  13,  the  deceleration  dropped  sharply 
as  the  penetrator  entered  Layer  3,  reaching  a level  of 
45-50  g's.  Again  assuming  that  this  level  would  remain 
fairly  constant  within  Layer  3,  and  using  a value  of  47 
g's,  it  is  estimated  that  the  penetrator  reaches  Layer  4, 
or  a depth  of  24  ft,  at  59.5  msec  and  with  a velocity  of 
101  m/sec,  or  330  ft/sec.  The  leapfrogging  approach  could 
have  been  used  to  start  a Phase  III  calculation  (covering 
penetrator  entry  into  Layer  4),  but  time  constraints  pre- 
cluded further  Integration  of  the  problem. 

Figure  14  shows  the  penetrator  velocity  vs  depth. 
Figure  15  shows  the  depth  of  penetration  vs  time.  Figures 
16  and  17  show  the  axial  force,  deceleration,  and  velocity 
vs  time. 


33 
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3 


4.2  PKNETRATOR-SOIL  INTERACTIONS  AND  SOIL  REPONSE 


4.2.1  Field  Plots 

Appendix  B contains  a series  of  field  plots  of  the 
computational  grid,  particle  velocities,  and  principal 
stresses  in  the  soil  target  medium  around  the  penetrator. 
These  plots  are  for  penetrator  deptli  j of  10  cm,  40  cm,  200 
cm,  280  cm,  and  528  cm. 


4.2.2  Friction 

In  addition  to  the  total  axial  force  on  the  penetra- 
tor and  its  deceleration,  Figure  13  shows  the  axial  force 
and  deceleration  due  to  friction  on  the  penetrator-soil 
interface.  Figure  18  shows  the  frictional  force  as  a per- 
centage of  the  total  axial  force  on  the  penetrator.  It  is 
seen  that  the  friction  rule  used  in  this  calculation  was 
responsible  for  a 40-60  percent  of  the  penetrator  decelera- 
tion in  all  the  layers  considered. 

Figure  19  shows  the  loading  distribution  along  the 
penetrator  nose  in  terms  of  the  normal  and  tangential 
(frictional)  stresses.  This  plot  is  made  with  the  nose  at 
200  cm  (6.6  ft),  entirely  within  Layer  1.  (Loading  distri- 
butions for  other  depths  are  given  in  Appendix  C.)  The 
normal  stress,  , falls  sharply  with  increasing  distance 
from  the  nose,  dropping  to  zero  at  about  38  cm,  or  15  in. 
However,  frictional  stress,  t , remains  fairly  level  along 
the  nose.  This  is  because  t is  limited  in  the  friction 
rule  by  the  yield  surface.  The  yield  surface  for  Layer  1 
is  fairly  flat,  as  seen  in  Figure  7. 
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LAyar  2 


Depth  of  Penetration I ft 

Figure  18.  Ratio  of  Frictional  Force  to  Total  Axial  Force  on 
Fanetrator  va  Depth  of  Penetration 


When  the  normal  stress,  , drops  below  zero,  the 
grid  can  open  to  represent  separation.  This  occurs  as 
soil  material  flows  around  tho  ogival  nose,  nearing  the 
tangoncy  point.  Tangency  with  the  cylindrical  body  is  at 

about  18.5  in.  In  Layers  1 and  2,  became  tensile  and 

target  material  separated  from  the  penetrator  at  about  38 

cm,  or  15  in.,  as  seen  in  the  plot  in  Figure  20.  After 

separation,  soil  may  subsequently  rebound  against  the  aft 
surface  of  the  penetrator.  This  did  m occur  in  the  cur- 
rent calculation  (even  if  the  deleted  cells  are  taken  into 
consideration),  but  it  should  be  expected  as  the  velocity 
drops  further.  While  normal  stresses  acting  against  the 
body  as  a result  of  such  rebound  will  presumably  be  small, 
they  may  be  applied  over  a large  area,  and  can  hence  lead 
to  significant  frictional  drag  at  late  times  in  the 
penetration. 

While  frictional  stresses  are  fairly  level  along  the 
penetrator  nose,  the  frictional  forco  distribution  is  much 
different.  Figure  21  shows  the  frictional  force  per  unit 
length  along  the  penetrator  in  Layers  1,  2,  and  3.  In  all 
layers,  friction  forco  peaks  several  inches  back  from  the 
nose.  The  frictional  force  is  still  quite  high  near  the 
point  where  the  target  soil  separates  from  the  penetrator, 
despite  the  fact  that  normal  stresses  drop  sharply.  This 
is  due  to  the  increasing  surface  area  of  the  nose  as  dis- 
tance from  the  tip  increases,  and  to  the  yield  surface  li- 
mitation on  frictional  stresses  in  the  region  where  normal 
stresses  are  high  (i.e,,  nearer  the  tip). 
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4.2.3  Prosauros  and  Stresses  in  Target  near  Penetrator 
Wose 

Figure  22  contains  vertical  profiles  of  pressure  in 
the  target  soil  along  the  axis  and  next  to  the  penetrator 
nose.  (The  plotted  pressures  are  from  the  first  row  of 
computational  cells,  which  are  attached  to  the  axis  ahead 
of  the  penetrator,  and  are  in  contact  with  the  penetrator 
surface  between  the  nose  tip  and  the  separation  point.) 
Three  profiles  are  given,  corresponding  to  representative 
depths  in  Layers  1,  2,  and  3.  Pressures  build  up  rela* 
tlvely  slowly  ahead  of  the  penetrator,  reaching  only  100* 
200  psi  10  cm  ahead  of  the  penetrator.  Peak  values  ('’1000 
psi  in  Layers  1 and  2)  occur  near  the  nose  tip.  There  are 
large  oscillations  as  unloading  starts  (probably  due  to 
the  finite  nature  of  the  grid  and  the  sharply  hysterotic 
unloading  model),  but  unloading  is  otherwise  fairly  linear 
to  the  separation  point. 

Plots  of  radial  stress  profiles  in  the  soil  outward 
from  the  penetrator  are  given  in  Figures  23  and  24  at 
penetration  depths  of  127  cm  (Layer  1)  and  300  cm  (Layer 
2).  This  type  of  plot  is  useful  in  assessing  the  adequacy 
of  the  zoning  employed.  The  sharp  stress  gradient  near 
the  nose  tip  Indicates  that  finer  zoning  in  the  radial 
direction  may  be  needed.  The  radial  zoning  employed  is 
probably  adequate,  however,  at  the  halfway  point  and  be- 
yond. A study  to  determine  the  zoning  requirements  for 
earth  penetrator  problems  is  planned  as  part  of  future 
work. 
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Figure  22.  Pressore  vs  Vertical  Distance  £m  Hose  Tip 
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Figure  2U.  Kedlel  Streii  ve  Redial  Dlatanoe  from 
Fenetratori  t - 21.0  maeo  (Layer  2) 
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4.3  COMPARISONS  WITH  EXPERIMENTS 


Preliminary  records  are  available  from  the  field  tests. 
The  records  which  are  shown  here  for  comparison  with  the 
calculated  results  are  smoothed  versions  of  the  uncorrected 
Sandia  traces,’  and  may  differ  somewhat  from  the  final  fil- 
tered and  corrected  curves  to  be  published  by  Sandia.' 

Figures  2S,  26,  and  27  compare  the  measured  and  compu- 
ted penetrator  deceleration,  velocity,  and  depth  histories. 
The  agreement  is  good.  The  calculated  decelerations  are 
somewhat  high,  but  the  ordering  of  decelerations  in  the 
three  layers  is  correct. 

There  are  a number  of  factors  which  may  have  contri- 
buted to  the  observed  differences: 

a.  Actual  site  profile  vs  Idealized  site  profile. 

The  experimental  accelerations  show  large  vari- 
ations within  each  of  the  layers  specified  in 
the  idealized  site  profile.  This  suggests  sig- 
nificant property  changes  within  the  layers. 

Also,  a preliminary  post-test  assessment  of  the 
site  has  Indicated  that  the  moisture  content  in 
the  top  4 ft  was  higher  than  had  been  assumed  for 
the  idealized  site  profile,  and  that  this  near- 
surface material  consisted  of  poorly-compacted  back- 
fill rather  than  undisturbed  soil  Cas  had  been  as- 
sumed in  the  constitutive  properties  recommended  by 
WES).^’  This  would  indicate  a softer  net.r- surface 
layer  and  presumably  lower  accelerations  in  that 
region.  Use  of  softer  properties  in  the  calcula- 
tions would  have  brought  the  calculated  values  into 
better  agreement  with  the  experimental  data  for  the 
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Ison  of  Conpated  and  Experiaeiital  Tine  Histories  of  Penetrator  Velocity 


Counted  awl  EjQteriaeital  Tiae  HistoEles  of  S^di  of  Penetration 


0-4  ft  depth  range.  Detailed  information  regard- 
ing the  actual  site  profile  will  be  forthcoming 
from  an  investigation  being  carried  out  by  WES. 


b.  Friction  rule.  Friction  played  a dominant  role 
in  the  numerical  solution,  as  has  been  seen  in 
I'igure  2.  Changing  the  coefficient  in  the 

friction  ruJ-*  from  .6  to  .4,  for  example,  would 
have  reduced  the  deceleration  level  in  Layer  1 
from  110  g's  to  89  g's.  Since  the  solution  re- 
sults (and  presumably  the  physical  phenomena) 
are  so  sensitive  to  the  friction  rule,  experi- 
ments and  more  careful  modeling  are  needed  to 
characterize  friction  under  the  loading  condi- 
tions involved  in  high  velocity  penetration. 

c.  Instrumentation  uncertainties.  The  output  from 
the  penetratOT  Instrumentation  may  be  subject 
to  errors  or  uncertainties  from  battery  excur- 
sions, inadequate  gauge  response,  and  body  vi- 
brations of  the  non-rigid  projectile.  The 
extent  of  these  problems  is  being  assessed  by 
the  experimental  investigators. 

A comparison  of  the  calculated  peak  pressure  vs  depth 
at  a radius  of  .45  m (1,48  ft)  from  the  test  centerline 
with  come  stress  data  obtained  in  the  test  is  shown  in 
Figure  28.  The  peak  pressure  vs  radius  in  Layer  2 is  com- 
pared in  Figure  29.  Unfortunately,  there  were  significant 
differences  between  the  calculational  conditions  (penetra.- 
tor  diameter,  mass,  and  velocity)  and  the  actual  test  con- 
ditions for  Sandia  tests  2 and  4 (as  indicated  in  the 
tables  on  these  figures),  so  direct  comparisons  cannot  be 
made.  However,  the  general  trends  appear  satisfactory. 


SECTION  V 


PROBLEM  AREAS 

In  carrying  out  the  finite  difference  rigid-body  cal- 
culation of  the  ONA  penetrator  problem,  a number  of  prob- 
lem areas  were  Identified.  Where  so  stated  in  the  follow- 
ing, these  areas  are  being  addressed  in  connection  with 
continuing  studies  of  earth  penetrators. 

a.  Large  distortions  occur  near  the  penetrator- target 
medium  interface,  due  to  radial  flow  around  the  advancing 
nose  and  to  shear  (friction)  along  the  sliding  surface. 
These  physical  distortions  produce  corresponding  distor- 
tions of  the  Lagranglan  computational  gr-ld,  necessitating 
periodic  realignment  of  the  grid  (l.e.,  local  rezonlng). 

If  done  carefully,  such  realignments  do  not  significantly 
degrade  a sblutlon,  but  they  are  time  consuming.  Automatic 
routines  are  therefore  being  Implemented  to  efficiently 
accomplish  periodic  realignment  of  the  distorted  grid  near 
the  penetrator.  (Grid  distortion  does  not  occur  in  a 
fixed  Bulerlan  grid,  but  a separate  set  of  problems  is  In- 
troduced due  to  radial  diffusion  of  target  material  pro- 
perties as  damaged  target  material  flows  through  the  fixed 
grid  out  and  around  the  penetrator  nose.) 

b.  Stress  gradients  are  affected  by  the  fineness  of 
the  computational  zoning.  In  order  to  determine  with 
reasonable  accuracy  the  stresses  in  the  target  medium  next 
to  the  penetrator  (where  gradients  are  steep),  relatively 
fine  numerical  zoning  may  be  required.  Studies  to  examine 
the  effects  of  zoning  on  penetrator  predictions,  and  to 
establish  zoning  criteria  for  reasonable  accuracy,  are 
being  undertaken. 
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c.  Friction  played  a very  important  role  in  this  DNA 
ponetrator  calculation,  accounting  for  at  least  half  of 
the  total  deceleration.  This  indicates  the  necessity  for 
experiments  and  more  careful  modeling  in  order  to  estab- 
lish realistic  friction  relationships  for  conditions  where 
metal  surfaces  slide  at  Siu''  velocity  and  under  high 
stresses  through  earth! <tu.  '..i  Such  experiments  are  being 
undertaken  under  other  DNA  .p’^ograms.  and  improved  friction 
modeling  is  being  undertaken  in  future  calculatlonal 
studies. 

d.  Information  regarding  the  properties  of  target 
sites  is  difficult  to  obtaih.  It  is  reasonable  to  assume 
that  this  will  always  be  the  case.  And  even  where  exten- 
sive subsurface  surveys  are  possible,  the  property  data 
are  uncertain  due  to  local  variations  and  to  in  aitu 
effects.  This  problem  is  quite  troublesome  if  one  is  at- 
tempting to  validate  a numerical  method  by  comparison  with 
a specific  test.  However,  exact  knowledge  of  target  pro- 
perties is  not  necessary  for  most  of  the  useful  applica- 
tions of  the  code;  l.e.,  for  studying  phenomenology  of 
penetration  in  different  basic  types  of  media,  or  for 
examining  the  effects  of  target  property,  penetrator  de- 
sign, and  Impact  parameters  upon  penetration  processes,  or 
for  establishing  the  design  environment  (stresses  and  ac- 
celerations] within  an  earth  penetrator  for  representative 
types  of  targets  and  Impact  conditions. 

e.  Computer  running  times  for  complete  penetration 
solutions  are  relatively  long.  Improvements  are  being  made 
which  will  substantially  reduce  these  times.  Nonetheless, 
code  solutions  will  probably  be  primarily  useful  for 
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analyzing  those  portions  of  penetration  events  where  load- 
ing on  the  penetrator  is  most  severe  or  is  changing  rapid 
ly  (e.g.,  while  the  penetrator  nose  is  burying  at  the  sur- 
face, or  while  the  penetrator  is  passing  between  different 
geologic  layers,  or  perhaps  just  before  it  comes  to  rest). 


SECTION  VI 

INTERNAL  STRESS  RESPONSE  OF  PENETRATOR 


In  a deep  penetration)  target  material  undergoes  large 
distortions  and  failure,  while  a successful  penetrator  de> 
sign  suffers  only  elastic  strains.  (A  sacrificial  nose 
section  might  be  an  exception).  In  penetration  analyses, 
the  projectile  can  therefore  be  assumed  to  be  a rigid  body, 
since  projectile  deformations  are  so  small  that  they  will 
not  significantly  interact  with  the  target  response  or 
penetration  dynamics. 

Rational  design  of  the  penetrator  interior,  however, 
requires  knowledge  of  the  design  environment,  specifically 
the  accelerations,  stress  histories,  and  shock  spectra  ex- 
perienced at  several  relevant  stations  within  the  penetra- 
tor. Numerical  solutions  in  which  the  penetrator  is 
treated  as  a deformable  body  provide  a means  for  obtaining 
such  Information. 

This  type  of  solution  is  Illustrated  in  Figures  30 
and  31,  performed  by  CRT  under  another  program.’  These 
figures  show  the  deformation  and  stress  field  occurring 
after  20  in.  penetration  into  soil  by  a 600  lb  projectile 
fired  at  2000  ft/sec. 

Finite  difference  solutions  of  this  type,  in  which 
the  penetrator  deformation  and  the  target  deformation  are 
treated  tog«th0v,  are  very  time  consuming  due  to  the  mar- 
kedly different  character  of  the  penetrator  and  target  re- 
sponses. Good  spatial  resolution  is  needed  in  the  penetra- 
tor in  order  to  resolve  structural  details,  and  the  wave 
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Flgura  31.  Principal  Straaa  Flald  with  Panetrator  at  20  In.  Depth  from 
Deformable  Body  Solution 


speeds  in  typical  penetrator  materials  are  relatively  high. 
These  factors  dictate  small  integration  time  steps.  On  the 
other  hand,  coarser  spatial  resolution  and  lower  wave 
speeds  in  the  target  media  would  permit  much  longer  time 
stops. 

Because  of  these  large  differences  in  the  integration 
time  step  and  in  the  magnitudes  of  the  distortions  involved, 
it  is  efficient  to  perform  separate  but  coupled  analyses 
of  target  penetration  and  of  penetrator  response.  Thus,  a 
rigid-body  analysis  is  first  performed  of  the  penetration 
in  the  target  media.  From  this  rigid-body  solution,  the 
penetration  history  and  force  loading  history  Imposed  on 
the  projectile  surface  is  determined.  This  information  is 
then  used  to  drive  a deformable  body  analysis  of  the  stress 
wave  response  in  the  penetrator  itself. 
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APPENDIX  A 


MATERIAL  MODELS 

Material  models  for  each  of  the  four  site  layers  were 

developed,  based  on  the  set  of  recommended  soil  properties 
furnished  by  WES.  Hysteretlc  elastic- ideally  plastic  models 
using  a Prandtl-Reuss  (non-aasociated)  flow  rule  were  employed. 
For  loading,  the  bulk  and  shear  moduli  were  made  functions  of 
the  volumetric  strain  (p) . For  unloading,  the  bulk  modulus 
was  made  a function  of  or  the  maximum  pressure 

and  the  shear  modulus  was  constant,  except  for  Layer  4, 
where  it  varied  with  the  unloading  bulk  modulus.  The  failure 
surfaces  were  fit  by  making  the  yield  condition  a function  of 
pressure. 

A.l  LAYER  1 

A. 1.1  Loading  Hvdrostat 

The  following  equations  were  doducini  for  the  loading  hy- 

drostat: 

P • a,y  U<0 


0<p<Pi 


P ■ P,  + y(u-P,)  + 


RFUp’  ■ F 


Ui<UiUc 
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This  model  was  extended  to  high  pressures  and  densities 
(beyond  that  furnished  in  the  representation)  by  assuming 
that  the  material  locked-up,  or  reached  void  closure,  at  a 
stress  of  '12  ksi  aiul  a strain  of  “38%  (Vj^  was  given  as  '37%) 
Beyond  void  closure,  a wcll-bohaved  asymptotic  form  is  em- 
ployed. The  values  of  the  constants  are; 


7.8933  X 10  ^ Mb 
‘.019894  Mb 
.38493  Mb 
-4.0302  Mb 
22.519  Mb 
-62.409  Mb 
67.438  Mb 
1.S3904  X 10‘®  Mb 
.38941 


S.2S46  X 10' 


.018319  Mb 

.689465  Mb 

1.0191  X 10“^  Mb 

7.74858  X 10'^  Mb 

i.4897  gm/cm^  (93  Ib/ft®) 


A . 1 . 2 Shear  Modulus  - Loadlni 


The  shear  modulus  function  for  loading  was  approximately 
determined  from  a graphical  assessment  of  the  variation  of 
Poisson's  ratio  along  the  uniaxial  strain  stress  path  (.stress 
difference  vs  pressure)  and  the  bulk  modulus  (as  derived  from 
the  hydrostat  equations).  Parametric  Integration  runs  were 
then  made  to  find  values  of  the  constants  giving  a good  fit 
to  the  stress  path.  A maximum  value  equal  to  the  unloading 
shear  modulus  was  imposed  to  prevent  energy  generation. 


3B(l-2vl 

TTl'+v)' 


^unload 


H tJf-i  y.  <■  I.. 


whore 


B - (y  + 1)  W 


and  V - .5  ■ .25  exp (- 6. 5u * 

Initial  values  o£  Young's  modulus,  the  bulk  modulus,  the 
shear  modulus,  the  constrained  modulus,  Poisson's  ratio,  and 
the  dilatatlonal  sound  speed  from  the  model  are; 

Eq'  « 17,173  psi 

Bjj  - 11,448  psl 

Gq  ■ 6,869  psi 

Mg  - 20,607  psi 

Cg  “ 1 , 013  ft/sec 
A . 1 , 3 Unloading  Model 

In  the  hysteretlc  region,  u<M^  , a model  giving  a fan 
of  constant  slope  P-u  unloading  paths  was  used  (slopes  are 
function  of  j 


■’m.x  * 


'^raax^'^c 


P « 


Pmaxi"c 


In  this  region,  the  unloading  bulk  modulus  varied  from  a mini- 
mum of  41,7  ksi  to  a maximum  of  430  ksi,  equal  to  the  bulk 
modulus  for  loading  at  void  closure.  The  unloading  bulk  mo- 
dulus was  made  a function  of  the  maximum  pressure, 
PmgjjCUmax^*  ^ constant  value  of  2S  ksi  for  the  unloading 
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shear  modulus  was  assumed,  The  initial  Poisson's  ratio  for 
unloading  was  «ot  equal  to  that  for  loading  (v^  ■ ,25),  giv- 
ing an  initial  alue  of  41.7  ksi  for  the  unloading  bulk  mod- 
ulus. K is  a function  given  by 


[exp  (bP, 
^max'*’^ 


max 


'^max  ^ 


c 


where  P^iax  obtained  from  the  loading  hydrostat  equations 
evaluated  at  . The  constants  are; 

Bo  “ 2.87277  x lO’^  Mb  b - 1000 

Bj  " .0229  Mb  ' 


For  , unloading  is  assumed  to  bo  reversible  (fol- 

lows the  load  path]) . 

Hydrostatic  tension  is  limited  by  imposing  a minimum 
value  of  pressure:  P,^in  " ‘12,26  psi. 

A . 1 . 4 Failure  Surface 

A yield  function  Y(P)  was  fit  to  the  failure  envelope 
suggested,  where  Y ■ at  failure.  Beyond  the  data  re- 
gion, an  asymptotic  form  was  used  to  smoothly  extend  the 
failure  surface  to  a limit  of  Y ■ 7S0  psi  at  infinity.  The 
functions  de.-ived  are: 

Y - 20  + 1.6313  P -12.26  <p<Q 

Y - 20  + 408. 14[exp(2. 4994  x lO'^’p) 

- exp  (-  3.747  x lO'^P)]  0<P<1420 
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Y - 600  + 150  p^:7s¥.°29 


P>1420 


Here  the  units  for  Y and  P are  psi. 

A . 1 . 5 Compaialive  Results 

Plots  of  stress  vs  strain,  stress  difference  vs  pressure, 
and  pressure  vs  volumetric  strain  for  uniaxial  strain  load- 
unlodd  paths  computed  with  the  model,  along  with  the  WES 
curves,  are  shown  in  Figures  A-1  to  A-6.  For  each  type  plot, 
the  first  figure  is  for  the  low  stress  region  and  the  second 
figure  is  for  the  entire  range  up  to  and  somewhat  beyond  void 
closure. 

A.  2 LAYER  2 

A. 2.1  Loading  Hydrostat 

The  following  equations  were  deduced  for  the  loading 
hydrostat j 


P • ajp 

p ■ 

i-1 


P ■ P + K (u-u  1 - (K  -K  lu* 
c *'max''^  '■  max 


1-exp 


Uc'P 


li^o 

0<y<Uj 


The  equation  for  the  fit  in  the  data  region  was  used  to 
extrapolate  the  hydrostat  up  to  an  approximate  void  closure 
point  corresponding  to  a stress  of  '48  ksi  and  a strain  of 
-41%  (V.  was  given  as  •'4L'‘.^.  Beyond  void  closure,  a well- 
behaved  asymptotic  form  is  employed.  The  values  of  the 
constants  are: 
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Strain  Stress  Paths  and  Failore  Sarface  for  Layer  1 (P  to  400  psl) 


Uj  - 2.27  X 10'^  Mb 
02  ■ *7.6853  X lO'^  Mb 
Oj  - .011947  Mb 
- - .034525  Mb 
Og  " .04S002  Mb 


Me  “ -7 


.0233697  Mb 
.689465  Mb 
3.15417  X 10"^  Mb 
1.4256  gm/cm^(89  Ib/ft^) 


A. 2. 2 Shear  Modulus  - Loading 

The  shear  modulus  In  loading  was  determined  by  a method 
similar  to  that  for  Layer  1. 

_ u4-r3Ba-2v)  ^ 1 


liiLl  r 
♦V)  * “unload 

m 


where 


B - (y+1) 


V « .5  - .3  exp(-.1687y*^^®) 


Initial  values  of  Young's  modulus,  the  bulk  modulus,  the 
shear  modulus,  the  constrained  modulus,  Poisson's  ratio,  and 
the  dilatational  sound  speed  from  the  model  are: 


Eg  > 5,926  psi 
Bg  ■ 3,292  psi 
G ■ 2,469  psi 


6,585  psi 


585  ft/sec 


A . 2 . 3 Unloading  Model 


In  the  hysterotic  rogion,  , a model  giving  a fan  of 

constant  slope  P-p  unloading  paths  was  used  (slopes  are  a 
function  of  ) J 


P ■ Vx  • '''’‘-“max) 


P • 


u <u 
^max  ^ 


c 

■c 


In  this  region,  the  unloading  bulk  modulus  varies  from  a min* 
Imum  of  74.3  ksl  to  a maximum  of  S76  ksl,  equal  to  the  bulk 
modulus  for  loading  at  void  closure.  A constant  represonta* 
tive  value  of  67  ksl  for  the  unloading  shear  modulus  was  as- 
sumed. This  value  was  determined  from  the  given  unload  path 
from  Oj  ■ 1000  pal  and  an  estimate  of  the  unloading  Poisson's 
ratio  at  that  point.  K is  a function  given  by: 


Bq+B,  [exp (bp 
^ » 


p TT 
*^max 


m.x)-‘l 


The  Constanta,  as  determined  from  the  given  unload  paths, 
are: 

Bq  ■ 5.12073  X 10’^  Mb 
Bj  ■ 7.8964  X 10'^  Mb 
b - 5.4326 

The  initial  Poisson's  ratio  for  unloading  is  .1532  from 
this  fit. 

For  P>Ug  , unloading  is  assumed  to  be  reversible  (follows 
the  load  path) . 


Hydrostatic  tension  is  limited  by  imposing  a minimum 
value  of  pressure:  - -12  psi. 

A . 2 . 4 Failure  Surface 

A yield  function  Y(P)  was  fit  to  the  failure  envelope 
suggested,  whore  Y ■ at  failure.  An  asymptotic  form 

was  used  to  smoothly  extend  the  failure  surface  in  the  high 
pressure  region.  The  functions  derived  are; 

Y ■ 16  ♦ jp 

Y - 416  + yCP-300) 

‘ 7. 7552  X 10“^  (P-300)^*®^^2 

Y « 1258.3  + 1647.87  P>1200 

Here  the  units  for  Y and  P are  psi. 

A . 2 . 5 Comparative  Results 

Plots  of  stress  vs  strain,  stress  difference  vs  pressure, 
and  pressure  vs  volumetric  strain  for  uniaxial  strain  load- 
unload  paths  computed  with  the  model,  along  with  the  WBS 
curves,  are  shown  in  Figures  A-7  to  A-12.  For  each  type  plot, 
the  first  figure  is  for  the  low  stress  region  and  the  second 
figure  is  for  the  entire  range  up  to  and  somewhat  beyond  void 
closure. 

A.  3 LAYER  3 

A. 3.1  Loading  Hydrostat 

The  following  equations  were  deduced  for  the  loading 
hydrostat : 


-12<P<300 


300<P<1200 
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Stress  Paths  and  Failure  Surface  for  I^yer  2 (P  to  400  psl) 


Volumetric  Strain  p 

rigure  A“12.  Hydrostatic  Load-Uhload  Paths  for  Layer  2 
(P  to  46  ksl) 


p 

U 

M<0 

p 

■ 

?!  + aj(y-u.)  + 

* 

CjCu-Ui)^ 

i - 1,2,5 

p 

u 

P,  - . A(u,-.0 

a 

p 

m 

Wc'Pl 

- t''«ax->'c’'‘* 

1-cxp 

m 

y>Uc 

The 

values  of  the  constants 

u re : 

“l 

m 

2.38661  X 10“'*  Mb 

- 0 

‘»1 

m 

-1.43803  X 10'^  Mb 

^2 

- .0547 

'^l 

m 

7.11694  X 10"^  Mb 

- .1013 

m 

2.86813  X 10’^  Mb 

W4 

- .1221 

m 

1.02409  X 10'^  Mb 

»^c 

- .14155 

= 2 

m 

6.39712  X 10'^  Mb 

■ .25 

*^3 

m 

7.99009  X 10'^  Mb 

•'c 

- .028085 

Mb 

to 

m 

9.96726  X 10‘^  Mb 

V 

max 

- .689465 

Mb 

^=3 

m 

•1.80756  X 10'^  Mb 

Pi 

- 0 

A 

m 

3.79269  X 10'2  Mb 

P2 

■ 1.37893 

X 10'®  Mb 

a 

m 

1.11227 

P3 

- 3.58522 

X 10’^  Mb 

Po 

m 

1.8S81  gm/cm^C116  Ib/ft^ 

) Pc 

- 1.27413 

X 10‘^  Mb 

A.  3. 2 Shear  Modtilus  - Loading 

The  shear  modulus  function  for  loading  was  developed 
from  an  analytic  fit  to  Poisson's  ratio  along  the  loading 
stress  path  and  the  bulk  modulus  as  derived  from  the 


hydrostat  equations.  A maximum  value  e(|ual  to  the  unloading 
shear  modulus  was  imposed  to  prevent  energy  generation 


“load  ■ 


3B(l*2v) 


p 

unload 


where 


B - Cu+1)  


and 


« dA 

„ . 

■ rrro 


where 


A ■ 0,-  o_ 
z r 


The  following  expressions  were  derived  for  the  derivative 


^ along  the  uniaxial  stress  loading  path: 


^ - S + 


0<,P<300  psi 


dA  , DB 
^ (E+P-300)^ 

where 

S - .66667 

B - “4.69643  x 10*^  psi'^ 
0 - Z. 07216 
D - 137.925  psi 
E « 688.145  psi 


P>300  psi 


Initial  values  of  Young's  modulus,  the  bulk  modulus, 
the  shear  modulus,  the  constrained  modulus,  Poisson's  ratio, 
and  the  dilatational  sound  speed  from  the  model  are: 


Ejj  - 3,115  psi 
- 3,462  psi 
Gq  ■ 1,154  psi 


M - 5,000  psi 
o ^ 


■ 447  ft/sec 

0 


A. 3. 3 Unloading  Model 

In  the  hysteretic  region,  u<:y^  , a model  giving  a fan 
of  constarit  slope  P-y  unloading  paths  was  used  (slopes  are 
a function  of  ) : 


'^max  * ^^^‘^max^ 


^max^^c 


P - Pc  + KgCy-y^) 


max-^c 


In  this  region,  the  unloading  bulk  modulus  varies  from  a min- 
imum of  35.8  ksi  to  a maximum  of  465  ksi,  equal  to  the  bulk 
modulus  for  loading  at  void  closure.  A constant  representa- 
tive value  of  IS  ksi  for  the  unloading  shear  modulus  was  as- 
sumed. This  value  was  determined  from  the  given  unload  path 
from  0^  ■ 1000  psi  and  an  estimate  of  the  unloading  Poisson's 
ratio  at  that  point.  K is  a function  given  by: 

. * Bi[.xp(b«„„)-l] 


U 

^max 


^max^^c 


< iVVWwir-.^ 


The  constants,  as  determined  from  the  given  unload  paths,  are; 

B ■ 2.47298  X 10‘^  Mb 
0 

- 1.02046  X 10*^  Mb 
b - 88.855 

The  Initial  Poisson's  ratio  for  unloading  is  .3165  from 
this  fit. 

For  unloading  Is  assumed  to  be  reversible  (follows 

the  load  path) . 

Hydrostatic  tension  Is  limited  by  imposing  a minimum 
value  of  pressure:  ■ ■ 29.5  psl  . 

A. 3. 4 Failure  Surface 

A yield  function  Y(P)  was  fit  to  the  failure  envelope 
suggested,  where  Y ■ failure.  An  asymptotic  form 

was  used  to  smoothly  extend  the  failure  surface  in  the  high* 
pressure  region.  The  functions  derived  are: 

Y - 22  + P *29.5<P<0 

Y ■ 22  + P - 4.0337  x o<P<300 

„ _ . 118.77(P-300) 

^ '6S'2\'J)7  * ^ 

Here  the  units  for  Y and  P are  psi. 

A . 3 . 5 Comparative  Results 

Plots  of  stress  vs  strain,  stress  difference  vs  pressure, 
and  pressure  vs  volumetric  strain  for  uniaxial  strain  load- 
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il.V  ulH.Ufc.  iu.' H^L«— 


j«sm: 


unload  paths  computed  with  the  model,  along  with  the  WES 
curves,  are  shown  in  Figures  A-13  to  A'16. 

A. 4 LAYER  4 


As  suggested  by  WES,  a model  using  a higher  value  for 
Poisson's  ratio  than  that  originally  recommended  was  formu- 
lated in  order  to  give  better  agreement  with  data  flrom  trl- 
axial  tests  with  small  constant  confining  pressures. 

A. 4.1  Loading  Hydrostat 

the  following  equations  were  deduced  for  the  loading 
'hydrostat: 


P ■ KqU  V<0 

P ■ P5  - ♦ a[«*p  0<u<Uj 


P - Pc  » 

The  values  of  the  constants  ar 

A - -3.572Z7  X 10*^  Mb 
a « -509.52 
a - 1.3 

Pg  - 1.43636  X 10"®  Mb 

A. 4. 2 Unloading  Model 

In  the  hysteretic  region, 
paths  were  made  a function  of 
the  value  of  attained: 


u>yc 

Kg  ■ 3.75514  X 10"^  Mb 

■ 1.02597  X 10’2  Mb 

Ug  • .0009 

pQ  ■ 1.9702  gm/cm^ 

(123  Ib/ft^) 

u<y_  , the  P-p  unloading 
the  hydrostatic  load  path  and 
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Stress  c , psl 


150 


§ 

T«d  *(^D 


s = s § 

I 

-*D)  BOUBjBjjfa  TBdTouTjj 
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load 


’’  ■ ’’load  • 


U >li 

*^max- c 


where  Pjoad^^^  loading  hydrostat  and  “ 3.7548 


X 10"^  Mb. 


For  > unloading  is  assumed  to  be  reversible 

(follows  the  load  path). 

Hydrostatic  tension  is  limited  by  imposing  a minimum 
value  of  pressure:  ■ -SO  psi  . 

■ 

A . 4 . 3 Shear  Modulus  - Loading  and  Unloading 

The  shear  modulus  for  both  loading  and  unloading  were 
made  variables » depending  on  the  bulk  modulus.  A con'^tant 
Poisson's  ratio  in  both  loading  and  unloading  was  employed. 


where 

end 


B - Cu+1) 


.497 


P(u)  is  the  appropriate  loading  or  unloading  hydrostatic 
equation. 

Initial  values  of  Young's  modulus,  the  bulk  modulus,  the 
shear  modulus,  the  constrained  modulus,  and  the  dilatatlonal 
sound  speed  from  the  model  are: 


E-  - 2,679  psi 


894.6  psi  c 


B, 


148.8  ksl 


Mg  » 150  ksl 


2,377  ft/sec 


''Hl1i»teiiiiii  

A . 4 . 4 Failure  Surface 

A yield  function  Y(P)  was  fit  to  the  failure  envelope 
suggested,  where  Y • failure.  The  functions  de- 

rived are: 

Y - 25  + ,5P  -S0<P<0  I 

Y - 50  + A(Pj  P)®  ♦ B(Pj-P)^ 

Y ■ 50 

where 

A ■ -2.58234  x 10'^®  psi 

B - -8.89951  X 10’^  psi 

a - 7.024 
- 150  psi 

Here  the  units  for  Y and  P are  psi. 

A. 4. 5 Comparative  Results 

Plots  of  atres.s  vs  strain,  stress  difference  vs  pressure, 
and  pressure  vs  volumetric  strain  for  uniaxial  strain  load- 
unload  paths  computed  with  the  model,  along  with  the  WES 
• curves,  are  shown  in  Figures  A-17  to  A-19. 


0<P<150 

P>1S0 


M'-i 
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N 200 


WES  Representation 


CRT  Model 


0 .02  .04  .06  .08  .10  .12 

Vertical  Strain  c^,  X (Engineering) 

Figure  A-17.  Uniaxial  Strain  Load-Unload  Paths  for  Layer  4 


Figure  A— 18.  Dni axial  Strain  Stress  Path  and  Failure  Surface  for  Layer 


APPENDIX  B 


DISTORTION,  VELOCITY,  AND 
PRINCIPAL  STRESS  FIELD  PLOTS 


A series  of  field  plots  of  the  computational  grid, 
particle  velocities,  and  principal  stresses  which  illus- 
trate the  soil  response  during  the  penetration  were  pre- 
pared. These  are  shown  for  the  following  penetration 
depths  and  times: 


Depth 

(m)  (ft) 


Figure  No.  of  plot  for 


Istortion 

(Grid)  Velocity  Stress 


.33 

1.32 

6.S6 

9.19 

9.83 

17.31 

17.58 


The  plots  of  the  grid  also  show  which  material  is  cur- 
rently failing;  i.e.,  is  on  the  yield  surface.  This  is  de- 
noted by  cells  containing  an  x or  -*•  , with  x indicat- 

ing a compressive  pressure,  P > 0 , and  + indicating 
hydrostatic  tension,  P < 0. 

The  velocity  vector  field  plots  show  the  direction 
and  magnitude  of  the  particle  velocity  at  each  lattice 


point  In  the  computing  grid.  (For  clarity  in  viewing  the 
soil  response,  the  velocity  vectors  of  the  projectile  are 
not  shown. j 


In  the  stress  field  plots,  the  principal  components 
of  the  stress  tensor  for  each  cell  are  shown,  as  follows: 
The  magnitude  of  the  two  principal  stresses  in  the  r-z 
plane  are  plotted  in  their  corresponding  principal  direc* 
tions.  The  third  principal  stre»s  (in  the  azimuthal  di- 
rection) is  plotted  along  the  line  bisecting  the  other 
two  principal  directions.  Vectors  pointing  to  the  right 
are  compressive,  to  the  left,  tensile.  An  example  of 
how  a stress  tensor  is  plotted  Is  sketched  below: 


Tensile 


Cell  Center 


Compressive 


Principal  Stresses 
in  r-z  Plane 


Principal  Stress 
in  0 Direction 
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Figure  B-4.  Grid  Configuration  with  Penetrator  at  40  cm  Depth 
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Figtire  B-8.  Particle  Velocity  Field  vith  Penetrator  at  f?00  cm  Depth 
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Figure  B-17.  Grid  Configuration  at  End  of  Solution,  with 
Penetrator  at  536  cm  Depth 
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APPENDIX  C 


FRICTION  STRESSES,  TARGET  DYNAMICS 
NEAR  PENETRATOR  NOSE,  AND  ENERGY  PARTITIONING 

Normal  and  tangential  (frictional)  stresses  applied 
to  the  penetrator  nose  surface  are  shown  in  Figures  C-1 
to  C-S  for  penetrator  depths  of  10,  40,  and  200  cm  (in 
Layer  1),  280  cm  (in  Layer  2),  and  S28  cm  (in  Layer  3). 
Generally,  these  plots  show  that  the  normal  stress  falls 
sharply  with  increasing  distance  from  the  nose,  The 
frictional  stress,  on  the  other  hand,  is  nearly  level  in 
Layers  1 and  3 (due  to  the  constraint  on  friction  imposed 
by  the  relatively  flat  yield  surface  in  these  layers).  The 
higher  yield  surface  in  Layer  2 allows  higher  frictional 
stresses. 

Velocity  components  of  the  target  soil  points  which 
are  in  contact  with  the  penetrator  nose  are  shown  in 
Figures  C*6  to  C*8  at  penetration  depths  of  200  cm  (Layer 
1),  300  cm  (Layer  2),  and  536  cm  (Layer  3).  In  these 
plots,  Vj^  and  are  the  velocity  components  normal 

and  tangential  to  the  penetrator  surface. 

Time  histories  of  the  kinetic  energy  of  the  penetra- 
tor and  the  kinetic  and  internal  energy  in  the  target  are 
shown  in  Figure  C-9. 
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Vertical  Distance  from  Moae  Tip,  cm 


Figure  C-: 


. Normed  and  Tangential  Stress  (fir) 
Distributions  along  Fenetrator  Nose 
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Figure  C-8>  Velocity  Components  of  Soil  In  Contact  with  Projectile  Nose 
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APPENDIX  D 


TIME  HISTORIES  OP  STRESS,  VELOCITY,  AND 
DISPLACEMENT  OF  TARGET  STATIONS 


Time  histories  oi'  stress,  velocity,  and  displacement 
were  obtained  at  a number  of  locations  in  the  soil  medium. 
The  stations  selected  for  the  code  output  correspond  to  grid 
points  nearest  to  the  layout  of  stations  suggested  by  WES. 
The  following  parameters  are  plotted: 

a.  Radial  stress 
Vertical  stress 
Hoop  stress 
Shear  stress 


r 

02 

"rz 


b. 


u 

e 

w 


u 

w 


Radial  particle  velocity 
Vertical  particle  velocity 

c.  Radial  displacement 
Vertical  displacement 

Positive  stresses  are  compressive.  Positive  u,  ii  is 
rudially  outward.  Positive  w,  dr  Is  vertically  down. 

The  locations  of  the  time-history  stations  in  Layer  1 
are  shown  in  the  following  sketch.  No  meaningful  time 
histories  were  recorded  in  Layer  2,  due  to  the  limited  pene- 
tration of  that  layer  by  actual  code  Integration.  Partial 
stress  and  velocity  time  histories  are  included  for  one  sta- 
tion (HZ)  in  Layer  3,  located  at  a radius  of  -.S  ft  and.  a 
depth  of  'I?. 5 ft. 

As  an  aid  in  interpreting  the  graphs,  the  time  of  arri- 
val of  the  projectile  tip  at  the  same  depth  as  the  station 
depth  is  indicated  on  the  graphs  by  a vertical  dashed  line 
labelled  t„  . 
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Figure  D-6 . Stress  Conponents  at  Station  3 
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Figure  0-11.  Stress  Components  at  Stetlon  A 
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Figure  D-12.  Velocity  Components  at  Station  4 
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